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Abstract 

Contact algorithm between different bodies plays an important role in solving collision problems. 
Usually it is not easy to be treated very well. Several ones for material point method were proposed 
by Bardenhangen, Brackbill, and SulskyQ Q, Hu and Chen[l^]. An improved one for three- 
dimensional material point method is presented in this paper. The improved algorithm emphasizes 
the energy conservation of the system and faithfully recovers opposite acting forces between contact- 
ing bodies. Contrasted to the one by Bardenhagen, both the normal and tangential contacting forces 
are more appropriately applied to the contacting bodies via the contacting nodes of the background 
mesh; Contrasted to the one by Hu and Chen, not only the tangential velocities but also the nor- 
mal ones are handled separately in respective individual mesh. This treatment ensures not only the 
contact/sliding/separation procedure but also the friction between contacting bodies are recovered. 
The presented contact algorithm is validated via numerical experiments including rolling simulation, 
impact of elastic spheres, impact of a Taylor bar and impact of plastic spheres. The numerical re- 
sults show that the multi-mesh material point method with the improved contact algorithm is more 
suitable for solving collision problems. 
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I. INTRODUCTION 



Phenomena with large deformation and/or large rotation are very common in nature, 
especially in fields of hypervelocity impact and explosion. Numerical simulations of such pro- 
cesses are necessary and challenging. The material point method(MPM) is an extension of 
FLIP 0, 3 which combines the strength of Eulerian and Lagrangian descriptions of the ma- 
terial, to the solid mechanics. The Lagrangian description is provided by discretizing each 
body by a collection of material points, and the Eulerian description is based on a background 
computational mesh. Information carried by the material points is projected on to the back- 
ground mesh where equations of motion are solved. The mesh solution is then used to update 
the material points. In Sulsky et aljsl, 0] a weak formulation of the MPM algorithm for solid 
mechanics is given and the method is framed in the terms of finite elements. The MPM 
combines the advantages of Eulerian and Lagrangian methods, which can avoid the distor- 
tion of Lagrangian mesh and track the boundaries of bodies. The method has been applied 
to the large strain problems 7, @], calculations with dynamical energy release ratefllj], frac- 
ture mechanics [12I ] , d ynamics failure^, 0], hypervelocity impact jsj, the thin membranes [lo| . 
granular materials [3, 0, 0, etc. 

In MPM, no slip contact between bodies is contained in the basic algorithm without 
additional cost. But at most cases separation or sliding may happen during the moving of 
bodies. A contact algorithm was presented by Bardenhagen, Brackbill and Sulsky to simulate 
the interactions of the grains of granular material 13]. In the algorithm, the contact may occur 
if the material points of different bodies are projected on to the same nodes of the background 
mesh, and the contact force is associated with the center-of-mass velocity. Bardenhagen's 
algorithm is linear in the number of grains and allows separation, sliding and rolling. With 
the contact algorithm MPM has been successful in simulating the large deformation of shear 
in granular material, having an advantage over traditional finite element methods(FEM) in 
that the use of regular grid eliminates the necessity to do costly search for contact surfaces. In 
order to apply MPM to stress propagation in the granular material, the contact algorithm is 
improved by Bardenhagen et al 14]. In Bardenhagen's improved contact algorithm, the normal 



traction is included in the contact logic to more appropriately determine the free separation 
criterion. 

To release the no-slip contact algorithm in MPM, a multi-mesh mapping scheme is pro- 
posed by Hu and Chen 18]. In the multi-mesh mapping scheme, each material lies in an 
individual background mesh rather than in the common background mesh. The meshing 
process of spur gears is simulated by Hu and Chen with their contact algorithm. To avoid 
interpenetration and allow separation in the gear meshing process, the normal velocity of 
any particle at the contact surface is calculated in the common background mesh, while the 
tangential velocity is found based on the corresponding information in respective individual 
mesh. With the proposed contact algorithm, Hu and Chen have successfully simulated the 
contact and separation of the gears. In their scheme, normal acceleration is set to be equal if 
particles of different bodies are mapped on the same node. But for some cases, the bodies may 
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separate although their particles are still mapped on the same nodes, which can cause energy 
dissipation during contact. In the contact algorithm by Hu and Chen, the friction between 
different bodies is completely ignored because the tangential velocities of different bodies are 
assumed to be independent. 

In this paper, an improved multi-mesh contact algorithm for three-dimensional MPM is 
proposed. In the present contact algorithm, the criterion of contact condition is similar to 
Bardengen's which ensures that the search for the contact of different bodies is fast, but 
the multi-mesh is used to calculate the normal and tangential velocities of different bodies. 
To avoid interpenetration the normal contact force is calculated at the contact surface and 
the Coulomb friction is applied in the tangential direction. Contrasting with Bardenhagen's 
algorithm, both the normal and tangential contacting force are more appropriately applied to 
the contacting bodies at the contacting nodes of the background mesh; Contrasting with Hu 
and Chen's contact algorithm, the normal velocities of different bodies are deal with separately 
just as the tangential velocities to not only ensure the contact/sliding/separation procedure 
can be simulated but also ensure the friction between different bodies can be applied. With 
the presented algorithm the total energy of the system is nearly constant during both elastic 
and non-elastic collision procedures, which shows numerical energy dissipation is little. 

This paper is organized as follows. The material point method is briefly reviewed in 
section [Til an d the new multi-mesh contact algorithm is illuminated in Section IIHI Several 
numerical examples are presented In section IIVI to validate the interaction between bodies 
with the contact algorithm. Numerical results obtained by the proposed contact algorithm 
presented are compared with those obtained by Bardenhagen's contact algorithm which show 
the proposed algorithm is more suitable in solving collision problems in which the numerical 
energy dissipation need to be very low. Section [V] concludes the paper with some remarks and 
observations. 



II. THE MATERIAL POINT METHOD 

The MPM is a particle method based on particle- in-cell (PIC) method in computational 
fluid mechanics. The method was initially developed for and has been successfully applied in 
problems involving large-deformation, large rotations of solid, etc. For continuum bodies, the 
conservation equation for mass is 

^+pV-v = 0. (I) 
And for pure mechanical problems the differential equation of balance is 

dv 

P^ = V-a + pb, (2) 

where p is the mass density, v is the velocity, a is the stress tensor and b is the body force. 

The formulation (j2j) is solved in a Lagrangian frame on a finite element mesh. The La- 
grangian formulation means that the momentum equation does not contain the convection 
term which can cause significant numerical error in pure Eulerian approaches. 
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In MPM, the continuum bodies are discretized with N p material particles. Each material 
particle carries the information of position, velocity, mass, density, stress, strain and all other 
internal state variables necessary for the constitutive model. Since the mass of each material 
particle is equal and fixed, Eq. is automatically satisfied. At each time step, the mass and 
velocities of the material particles are mapped onto the background computational mesh(grid). 
The mapped nodal velocity Vj is obtained through the following equation, 

Yl m 'J w i = Yl m pVp N i ( x p) (3) 

3 P 

where m p , v p and x p are the mass, velocity and position of particle p, separately. Ni is the 
element shape function, i and j indexes of node. For three-dimensional problem, a eight-node 
cell is employed with the shape functions given by 

^=i(l + £&)(l+Wi)(l + CCi) (4) 

where £, rj and £ are the natural coordinates of a material particle in the cell along the x-, y- 
and z-directions, rji and Q are the natural coordinates of the node i in the cell along the 
three directions. 

In the Eq. ([3]), the consistent mass matrix, my, is 

rriij = Y m p N i( x p) N j( x p) ( 5 ) 

p 

In practice, we generally replace with a lumped, diagonal mass matrix so that Eq. ([3l) 
becomes 

rriiVj = y^m p v p iVj(xp) (6) 

p 

where lumped mass is 

mi = y]m p iVj(xp) (7) 
p 

After the information is mapped from material particles to mesh nodes, the discrete for- 
mulation of Eq. (j2j) on the mesh nodes can be obtained, as described below. 

The weak form of Eq. (j2J) can be found, based on the standard procedure used in the 
finite element method, 

I p5w ■ dv/dtdfi + / 5(vV) • crdfi - f Sv-tdT- f p5v ■ bdQ = . (8) 
Jn Jq Jvt Jn 

where Q is the domain to be solved, Ft is the traction boundary, a is the stress tensor, t is 

the external force vector and b is the body force vector. 

Since the continuum bodies is described with the use of a finite set of material particles, 

the mass density can be written as, 

p( x ) = Y "vK x - x p) ( 9 ) 
p=l 
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where 5 is the Dirac delta function with dimension of the inverse of volume. The substitution 
of Eq. (j5J) into Eq. (??) converts the integral to the sums of quantities evaluated at the 
material particles, namely, 

= (f,) int + (f,) ext (10) 

where m ; is the lumped mass, (fj) mt and (fj) cxt are the external force vector and internal force 
vector which read separately 

N p 

(f 4 ) int = -^ m P a-(ViV 4 )/p p , (11) 
p 

N p 

(f 4 ) ext = Nib p + f? (12) 
P =i 

where the vector ff is the contact force which is the external nodal force not including the 
body force and is illustrated in the following section. 

An explicit time integrator is used to solve Eq. (jTUl) for the nodal accelerations, with the 
time step satisfying the stability condition. The critical time step is the ratio of the smallest 
cell size to the wave speed. After the equations of motion are solved on the cell nodes, the 
new nodal values of acceleration are used to update the velocity of the material particles. 
The strain increment for each material particle is determined with the use of gradient of the 
nodal basis function evaluated at the material particle position. The corresponding stress 
increment can be found from the constitutive model. The internal state variables can also be 
completely updated. The computational mesh may be discarded, and a new mesh is defined, 
if desired for the next time step. As a result, an effective computational mesh could be chosen 
for convenience. 



III. THE CONTACT ALGORITHM 

The MPM with a natural no-slip contact algorithm is based on a common background 
mesh. As a result, it is impossible to separate the contacting bodies. Bardenhagen et al. [3, [jjj 
have proposed a contact algorithm in which the contact between bodies is handled when the 
velocity field of an individual particle differs from the single, center-of-mass velocity field in 
the cell containing contacting particles. Their contact algorithm was incorporated into the 
MPM to simulate the interactions in granular materials based on the velocity field. 

In this section, we will improve the multi-mesh contact algorithm which recovers more 
faithfully the opposite acting forces between contacting bodies. In MPM, several bodies may 
be mapped on to the same nodes of the background mesh, so it is necessary to define a multi- 
value of velocity and mass on every node. In practice, it is impossible that the number of values 
defined at one node is as many as the number of bodies, otherwise the memory of computer 
will be too much wasted if there are thousands of bodies to be simulated. In this paper we 



define four values on every node. That is to say, there are at most four bodies mapped on 
to the same nodes, although there can be thousands of bodies in the whole domain. In that 
case, each node has a mesh mass mf and momentum Pf associated with it, where g ranges 
from one to four and the mesh velocity vf can be obtained from the mesh momentum and the 

mass, 

vf = Pf/mf (13) 

Note that if the mesh mass mf is close to zero, the obtained mesh velocity maybe singular 
which will cause error during the calculations. In this paper, to avoid the singularity, the 
shape function is altered, if £, rj or £ is small than —0.99 or larger than 0.99 that means the 
material point is too close to the node, £, r\ or ( is adjusted to —0.99 or 0.99. Since the shape 
functions have compact support, only those nodes in the vicinity of the bodies will have a 
meaningful velocity and the body velocity at other nodes will be zero. 

Obviously, if the momenta of two bodies are projected on to the same node, the contact 
may occur and the contact between bodies r and s is directed by comparing the fields vf and 
vf which are determined by using mass weighting given in Eq. ffTB"]) . 

(vf - vf) ■ nr > 0, (14) 

where nf s is the unit outward normal at node % along the boundary. Multiply Eq. (1T4"]) with 
mfmf, it can be written as, 

(mfPf - mfPf) ■ > (15) 

If Eq. ( !T4|) is satisfied, the velocities of body r and body s are adjusted to new values vf and 
vf so that 

vf • = vf • n[ s (16) 

holds. That is, the normal components of velocity of body r and body s are set to be equal. 
Eq. ({16]) can also be written as 

mfPf • n[ s = mfPf • n?. (17) 
As a reasonable contact algorithm, the momentum is required to be unaltered, i.e., 

(P[ + P- ) • n[ s = (P[ + P?) ■ nr (18) 
From Eqs. ( fTTl) and ([TBI the updated mesh momenta are obtained, 

Pf = Pf - (mfPf - mfPf) • <X7K + (19) 

Pf = Pf + «Pf - mfPf) • <X S /K + <)■ (20) 
So the updated mesh velocities are: 

v[ = v[ - mf (v[ - vf ) • nJXVK + (21) 

vf = vf + m[(v[ - vf) • <X7K + <)• (22) 



Especially, if body s is a rigid wall, we set the value of mf much larger than that of m£. Thus 
Eq. (J2TJ and Eq. ([22} can be reduced to: 

v[ = v[ - (vT - vj) ■ (23) 

vf = vf, (24) 

Obviously, the velocity of rigid body s is not altered during the contact. 

The equations fT2~Tj) and fT22l) determining the velocities are identical to that by 



Bardenhagen 13] and in practice they are same, but the calculation of the normal and tan- 
gential contact forces makes the difference which will be described as following are different. 

Once bodies r and s contact, they move together along the normal until they separate 
when the contact condition expressed in Eq. (fl5l) is not satisfied. So the acceleration along 
the normal of body r is equal to that of s during the course of the contact. That is 

a[ • <• = a* • nj- (25) 

where a£ and af are the accelerations of bodies r and s at node i, respectively. They can be 
obtained from the Newtonian second law, 

m[(a[ ■ nH = fT' int " < s - /T° r (26) 

m|(a| • nH = • n? + /T (27) 

where /™ or is the normal contact force between body r and body s at node i, which can be 
obtained from Eq. (|25l) (j^i 



fr = «ff nt - <r nt ) • <vk + o (28) 

Note that the normal contact force must be nonnegative. So once / 4 nor is negative, it is 
set to be zero. That is 

h "jo, *<o • (29) 

where \1/ = (m|fj r ' int — m[fj S ' mt ) • n^. For the cases /j nor in Eq. ( J28|) is positive, the contacting 
bodies at time t may still contact in the next time step although the criterion of contact in 
Eq.ffToT) is not satisfied, so the contact condition should be applied in the next time step. 
The contact force in the Bardenhagen's contact algorithm is calculated as 

= ml [(v< - v[) ■ n-] /At (30) 

where Vj is the center-of-mass velocity at node i. Actually, the normal contact force may still 
be very large even if the normal velocities of different bodies at contact nodes are same during 
the course of contacting. The normal contact force calculated by fl30|) is not physical and may 
cause numerical energy dissipation which will be shown in the next section. 

When without friction, the contact algorithm has been finished up to now. In the case 
with friction, the frictional slip is accomplished by adjusting the tangential component. To 
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apply Coulomb friction, we first calculate the force necessary to cause the bodies to stick 
together completely. Again, the comparison of the mesh velocity of body r to that of body s 
provides exactly the correct constraint for no-slip contact if body r and body s contact, the 
relative tangential velocity is 

(v[ - vf) ■ < (31) 
where s£ is the unit tangential at node i along the boundary, 

si = ((v[ - v!) - (vl - v|) ■ nrn[ s )/| ((v[ - v?) - « - yf) • n>r) I (32) 

To get an appropriate frictional force allowing slip, we start from the non-slip condition. 
The tangential velocities of body r and body s must be set to be equal after one time step 
At. Suppose the tangential contact force is /j tang . It should satisfy 

(P[ " + (ff" • Bf - fr s )At) jm\ = (Pf • af + (ff nt • s[ + /r s ) At) /< (33) 

Then the needed constraining tangential force / 4 tang is 

/f ns = (KP[ - mJP?) + Kf[' int - m[ff nt )At) • s[/(K + mJ)At) (34) 

The expected frictional force should equal to /* ang if the latter is small, and be proportional 
to the magnitude of the normal force and independent of the contact area if f i ans exceed 
a specified value. That is to say, the frictional force just balances the tangential force to 
prevent relative tangential motion when the latter is small. When the latter is large, we 
limit the frictional force to have a magnitude less than it to allow tangential slip between 
the contacting bodies. The direction of the frictional force is chosen as in ( |32l) to oppose the 
relative motion. Putting these requirements together yields, 

// ric = min(/i/r,/r S )- (35) 

where \i is the coefficient of friction. To complete the formulation of the contact algorithm, a 
value n[ s of the normal at node i of the computational mesh for the contacting bodies r and 
s is still needed. As an approximation, the following algorithm is presented to determine the 
normal value, 

1. If bodies r and s contact at node i, initialize vectors V r and V s as zero. 

2. Search within the eight cells(for three-dimensional cases) around the node. If cell j 
possesses particles belonging to body r (or s), calculate the difference of coordinates 
between the node i and the center of cell j, x, — x. 3 c , then add the difference to vector 
V r (or V s ). 

3. Calculate the difference between the vectors V r and V s , then set the difference as the 
value of n[ s . Finally, unitize n[ s . 

Finally, we summarize the material point method with the new multi-mesh contact algo- 
rithm presented in this paper as follows: 

R 



1. Get the particle mass m T , position x r , velocity v r , density p r , and stress a; form the 
lumped mass matrix(Eq. (j7|)) and nodal momentum(Eq. (jSJ)). 

2. Loop over the mesh nodes, if two bodies contact at node i, adjust the nodal momenta 
of the contacting bodies(Eq. (fl9l). Eq. (|20l)). 

3. Calculate the rate of the deformation gradient for each particle, compute the increment 
of strain using an appropriate strain measure and solve constitutive equations to update 
the stress, a p . 

4. Form the internal force (Eq. ffTTT) ) . Calculate the contact force between bodies and form 
the external force (Eq. (|12p ) 

5. Solve the momentum equations for the nodal accelerations and get the velocity in a 
Lagrangian frame: 

m t [ Vl \ t+At - Vi | t ] = At[(f,) int + (f,) cxt ] (36) 

6. Update the solution at the material point by mapping the nodal values using the element 
shape functions. Positions and velocities are updated according to 

x p | i+ At = Xp|t + At 2j Vi\t+AtNi(x p ) (37) 

i 

and 

Vp| t+ Ai = V p |t + ^2 t v il<+At - v i |t]JV i (x p ) (38) 

i 

7. Define a new finite element mesh if necessary, and return to step 1 to begin a new time 
step. 

IV. NUMERICAL SIMULATION 

Numerical simulations presented in this section are carried out in three dimension. The 
first set of simulation involves a cylinder rolling on an inclined rigid plane and is meant as 
simple illustration and validation of the friction algorithms presented by Bardenhagen et al. 
and by us. The second set involves the collision of two elastic spheres and is to examine the 
efficiency of the multi-mesh contact algorithm proposed in this paper. The third set involves 
a copper Taylor bar impacting to a rigid wall. The last example is to simulate the process of 
the collision between four identical spheres. The last two examples examine the conservation 
of energy during the collision is checked. 

A. Rolling simulation 

Fig Q] shows the plane geometry for a computation with a cylinder on an inclined plane. 
In this example, the plane inclined at an angel 9 with respect to the horizontal line, while the 
gravity g is vertically downward. 
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FIG. 1: Geometry for simulations of a cylinder on an inclined plane 



A rigid cylinder on an inclined surface will roll, or slip depending on the angle of inclination 
and friction coefficient. Specifically, if tan# > 3/i, the cylinder will roll and slip; Otherwise 
the cylinder will roll without slipping, where /i is the coefficient of friction. For an initially 
stationary, rigid cylinder, the ^-component and the center-of-mass position as a function of 
time, x cm (t), is given by 



In Eq. ( 139]) . xq is the a;— component of the initial center-of-mass position, and |g| is the 
magnitude of the gravitational acceleration. 

Simulation is performed with a cylinder that has the radius R = 40mm, thickness t = 
20mm, and gravitational acceleration with magnitude 10m/s 2 . The computational domain is 
cubic with side length length 700mm, 150mm and 40mm, respectively. The cell size is 10mm so 
there are only eight and two computational elements across the diameter and the thickness of 
the cylinder, respectively. The simulation involves a elastic, deformable cylinder with elastic 
modulus 1.24MPa, Poisson ratio 0.35 and density 8.0 x 10~ 3 g/mm 3 . The inclined plane is 
discretized as a rigid body with 70, 15 and 4 material points, respectively, and there is only 
one material point in one cell. 

Fig. [2] shows the center-of-mass position of the cylinder as a function of time for three 
values of angel of inclination, 9 = tt/12, 9 = tt/6 and 9 = 7r/4, respectively, and the coefficient 
of friction fixed at /x = 0.5. The symbols represent simulation results, and lines represent 
analytical ones. Fig. [2(a) shows the results of our contact algorithm while Fig. [2j(6) shows 
those of contact algorithm by Bardenhagen. For cases with large inclination angle the results 
of both contact algorithms agree well with analytical solutions. But when the inclination angle 
is small, results of our contact algorithm are much better. 

In the next test, the value of angle of inclination is fixed at 9 = tt/6 and the coefficient 
of friction is varied, \i = 0.1 and \i = 0.5. Fig. [3] shows the center-of-mass position of the 
cylinder as a function of time for each simulation and the corresponding exact solution for a 
rigid cylinder. The computed results for the deformable cylinders agree well with the analytical 
solutions, and as before, the computed curves obtained with our contact algorithm are much 
more closer to the analytical curves than those by Bardenhagen's algorithm. 

Fig. H] shows simulation results for different mesh sizes. In this test, the angle of inclination 
is fixed at 9 = ti/6 and the coefficient of friction is fixed at \x = 0.5. The side length of the 
cubic computational elements is varied, 40mm, 20mm and 10mm, respectively. Fig. H^a) 
shows the results of our contact algorithm while Fig. 2(6) shows those of Bardenhagen's. It is 
clear that the simulation results converge to the analytical ones with the decrease of the mesh 
size. The agreement between results by our contact algorithm agree better with analytical 
ones than those by Bardenhagen's scheme in the later time. 




xo + ||g|£ 2 (sin0 — //cos0), tan# > 3/i (slip), 
xq + ||g|t 2 sin 9, tan 9 < 3/i (stick), 
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FIG. 2: Center-of-mass position for deformable cylinder vs time. The coefficient of friction, fi = 0.5, 
the angles of inclination, 6 = tt/12, ir/6 and ir/4, respectively. The symbols represent simulation 
results while lines represent analytical ones, (a)our contact algorithm, (b) contact algorithm by Bar- 
denhagen. 




FIG. 3: Center-of-mass position for deformable cylinder as a function of time. The angle of in- 
clination, = 7r/6, and coefficient of friction, jjl = 0.1, 0.5, respectively. The symbols represent 
simulation results while lines represent analytical ones, (a) our contact algorithm, (b) algorithm by 
Bardenhagen. 

B. Impact of elastic spheres 

In this section, the impact of two elastic spheres is simulated to test the conservation of 
the energy during the impact with the contact algorithm. The variables are all dimensionless 
in this example. The spheres start from the left side and the center of the computational 
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FIG. 4: Center-of-mass position for deformable cylinder as a function of time. In the simulation 
angle of inclination is 9 = tt/6, coefficient of friction is /i = 0.5 and side lengths of cells are 40mm, 
20mm and 10mm, respectively, (a) our algorithm, (6) algorithm by Bardenhagen. 

FIG. 5: Snap-shots of impact of two elastic spheres obtained by our contact algorithm. From up to 
down, the corresponding times are 0ms, 37.8ms and 75.0ms, respectively. 

domain with initial velocities (0.1,0,0) and (0,0,0), respectively. The computational domain 
is a cube whose sides along the x, y and z direction are 40, 20 and 20, respectively, and cubic 
meshes are used with side length Ax = Ay = Az = 0.5. Eight material points are used per 
element. The spheres have a radius of 4, Young's modulus of 1000, a Poisson's ratio of 0.3 
and a density of 1.0. The distance between the center of first sphere and that of the second 
sphere is 14. The simulation is run up to a final time t = 80. 

The results from our contact algorithm are shown in Figj5j As a comparison, the results 
from the contact algorithm by Bardenhagen are shown in Figj6l In order to show the moving 
of spheres clearly, only the central layer of the 3D configuration is shown. The nonlinear 
constitution for large-deformation is used. In these figures, (a)show the two spheres at time 
t = 0.0 when they just begin to travel through the grid, (fr)show the impact of spheres at 
time t = 37.8 and (c)show the spheres at time t = 75.0. From FigJU (c) we find two spheres 
separate, and the left one is almost immobile and the right one moves with nearly the same 
kinetic energy as the initial kinetic energy of left one. But in Figj6] (c) the two spheres move 
forward together, which is unphysical. 

Figl7]shows the kinetic, potential and total energies as a function of time in which (a) shows 
the results of our contact algorithm while (b) shows those of Bardenhagen's. From Fig. [7(a) 

FIG. 6: Snap-shots of impact of two elastic spheres obtained by Bardenhagen's contact algorithm. 
From up to down, the corresponding times are 0ms, 37.8ms and 75.0ms, respectively. 
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FIG. 7: Energy evolution of elastic cylinder impact. The results with the contact algorithms presented 
in this paper and by Bardenhagen are shown in (a) and (6), respectively. 

we find the kinetic energy decreases during the impact and recovers after the spheres separate. 
The potential energy (broken line) begin to accumulate upon impact of spheres, reaches its 
maximum value at the point with maximum deformation during impact, then decreases to 
a small mount associated with the free vibration of the spheres after separation. The total 
energy(solid line) is approximately constant. In Fig. UXb), the kinetic energy decreases during 
the impact but doesn't recover to the original one; the potential energy begin to increase when 
the contact begins, then reaches its a maximum value, and does not decrease. Total energy 
shown in Fig. W(b) is not constant, which shows a strong numerical dissipation during the 
course of impact. 



C. Impact of a Taylor Bar 

The classical Taylor bar problem is considered. This is a commonly simulated problem 
and is often used as a benchmark for transient dynamics computer codes. A copper bar of 
radius R = 3.8mm and length L = 25.4mm impacts on a rigid, frictionless wall with an 
initial longitudinal velocity of 190m/s. The material is modeled as elastoplastic with Young's 
modulus E = 117GPa, Poisson v = 0.35, the yield stress is a y = 0.157MPa and linear 
hardening is assumed with H = 0.425MPa. The material has a density of po — 8.93g/cm 3 . In 
order to compare the computed results to those of experiments, we use the following estimation 



of error given by G. R. Johnson 19j: 
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A = o [ -T 1 + + ^77^ (40) 
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FIG. 8: The sketch figure of the impact of a Taylor bar 
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FIG. 9: Energy evolution during impact of Taylor bar 

where L and D are the length and diameter of the bottom after the impact, respectively, as 
shown in Fig. [HJ W is the diameters of the layer which is 0.2Lo to the bottom. 

The bar moves within the cubic domain [—10.4,10,4] x [—10.4,10.4] x [—2,36], meshed 
by 30 x 30 x 50 elements. For the initial construction of the bar, there are 8 material points 
in every cell, and for the rigid wall there is only one material point. The unit of coordinate 
is millimeter. The terminal time is 80/xs. Fig. [9] shows the kinetic, potential and total 
energy as a function of time, where the contact algorithm presented by us(Fig. [9](a)) and by 
Badenhagen(Fig. Mb)) are used. In Fig. E^a), the kinetic energy decreases during the impact 
and is totally converted to potential energy at the end. The total energy is constant during 
the whole time. In Fig. Myb), the energy is dissipated during the impact. 

Table [I] shows the comparison between the computed results and experimental ones, where 
MPM1 is MPM with our contact algorithm and MPM2 is MPM with contact algorithm 
presented by Bardenhagen. The results by MPM1 agree better with experimental ones. 

Fig. [10] shows the final particle configuration, colored by contour values of equivalent 
plastic strain obtained with MPM1. Fig. fTOT a) shows three-dimensional view while Fig. \W\ b) 
shows the center layer of Fig. [TOT a) vertically to the rigid wall. 



FIG. 10: The final particle configuration of the Taylor bar 
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TABLE I: The comparison between the computed results and experimental ones 





L(mm) 


D(mm) 


W(mm) 


A 


Experiment 


16.2 


13.5 


10.1 




MPM1 


16.15 


13.21 


9.63 


0.071 


MPM2 


16.25 


11.96 


9.42 


0.184 
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FIG. 11: Energy evolution during impact of four copper spheres. 
D. Impact of plastic spheres 

The last example simulates the impact of four identical copper spheres with the contact 
algorithm presented in this paper. The material parameters of spheres are the same as those 
of last example. The radius of the spheres is 10mm. Initially, one of the spheres locates 
at (0.0,0.0,25.0) and travels with a velocity — lOOm/s parallel to the z axis; The other three 
spheres locate at (10, -5.7735, 0.0), (0.0, 11.547, 0.0) and (-10.0, -5.7735, 0.0) are at rest. The 
unit of coordinate is millimeter. 

The spheres moves within the cubic domain [—50,50] x [—50,50] x [—50,50], meshed by 
50 x 50 x 50 elements. Fig. [11] shows the kinetic, potential and total energy as a function 
of time. The kinetic energy decreases during the impact and part of the kinetic energy is 
converted to potential energy during the impact. The total energy is a constant during the 
whole time. 

FigfTWaWc) show the particle configurations of different time, colored by contour value 
of equivalent plastic strain. From blue to red the plastic strain increases and correspondingly 
the possible temperature increment becomes larger. Fig. [TWa) shows the initial particle 
configuration at t = 0/is. Fig. [T2T 6) shows the particle configuration at t = 80/is when the 



1 r> 



FIG. 12: Snapshots of impact of four copper spheres. From black to white the plastic strain increases 
and correspondingly the possible temperature increment becomes larger. [a)t = O^ls; (b)t = 80/Us; 
and (c)t = 1ms 

upper sphere just impact to the lower three ones. Fig. [l~2f c) show the particle configuration 
at t = 1ms when they separate. 

V. CONCLUSION 

A new multi-mesh contact algorithm for three-dimensional material point method is pre- 
sented. The contact algorithm faithfully recovers the opposite acting forces between colliding 
bodies. Collision procedures between regular bodies and/or rigid bodies can be treated within 
the same framework. A multi-value of momentum and mass is defined on every node to de- 
scribe the contact/sliding/separation procedure. Both normal and tangential velocities of each 
particle at the contact surface are calculated in respective individual mesh. A Coulomb fric- 
tion is applied to describe the sliding or slipping between the contacting bodies. The efficiency 
of the contact algorithm is linearly related to the number of the contacting bodies because the 
overlapped nodes are labeled by sweeping the material particles of all bodies when the nodal 
momentum and mass are formed in every time step. 

Numerical simulation shows that our contact algorithm possesses high accuracy and low 
numerical energy dissipation, which is very important for solving collision problems. 
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